//==============================================================================
// Sequence.cc
// ----------------------------------------------------------------------------
// Class that represents a generic sequence and general methods for 
// working with them, including reading and writing sequences to 
// files.  
//------------------------------------------------------------------------------
// $Id: Sequence.cc,v 1.2 2000/02/01 23:42:47 kdurbin Exp kdurbin $
//------------------------------------------------------------------------------
// $Log: Sequence.cc,v $
// Revision 1.2  2000/02/01 23:42:47  kdurbin
// Added InitFromCString and ConcatenateCString
//   
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "Sequence.hh"

void
Sequence::ComputePadded2Depadded()
{
  mPadded2Depadded.resize(size());

  int depadIdx = 0;
  for(int i = 0;i < size();i++){
    char b = (*this)[i];
    if ((b != '*') && (b != '-')){
      depadIdx++;
    }
    mPadded2Depadded[i] = depadIdx;
  }
}


void
Sequence::ComputeDepadded2Padded()
{
  // Compute pad2depad first...
  if (mPadded2Depadded.size() != size()){
    ComputePadded2Depadded();
  }

  // Determine how big the depad array needs to be...
  int maxDepadIdx = mPadded2Depadded[size()-1];

  // Since it's one based, the size is the same as the maxIdx
  // However, the array will have an unused 0 slot, so still +1 for size.
  mDepadded2Padded.resize(maxDepadIdx+1);

  int lastDepad = -1;
  for(int i = 0;i < mPadded2Depadded.size();i++){
    int thisDepad = mPadded2Depadded[i];
    if (thisDepad != lastDepad){
      mDepadded2Padded[mPadded2Depadded[i]] = i; 
    }
    lastDepad = thisDepad;
  }
}

// Basically implements the following mapping:
//                   1         2
// -3 -2 -101234567890123456789012345 Padded
//         ACTG*GATCC*C*AAT*CATAC 
// -2  -1 012344567899001233456789012 Depadded
//                    1           2
// Note that coordinates can go beyond the range of the sequence itself!! 
// 
// 
int 
Sequence::Depadded2Padded(int DIdx){
  if (mDepadded2Padded.size() == 0){
    ComputeDepadded2Padded();
  }

  // Less than 1 because 1 is the leftmost sequence base. 
  // In what follows, we are looking at size() -1 because we have
  // an unused 0 location in the array.  So, our elements 
  // 1...10 are in an array that is 11 in size. 
  if (DIdx < 1){
    return(DIdx-1);
  }else if (DIdx < mDepadded2Padded.size()-1){
    return(mDepadded2Padded[DIdx]);
  }else{
    int lastDepad = mDepadded2Padded.size()-1;
    int diff = DIdx - lastDepad;
    return(mDepadded2Padded[lastDepad] + diff);
  }
}


// Basically implements the following mapping:
//                   1         2
// -3 -2 -101234567890123456789012345 Padded
//         ACTG*GATCC*C*AAT*CATAC 
// -2  -1 012344567899001233456789012 Depadded
//                    1           2
// Note that coordinates can go beyond the range of the sequence itself!! 
// 
// 
int  
Sequence::Padded2Depadded(int PIdx){
  if (mPadded2Depadded.size() != size()){
    ComputePadded2Depadded();
  }

  // Special case when padded index is less than zero.  Then we have this mapping:
  //
  // 
  //                        A G G G T A ...CAT
  //                                       11111
  // Padded        -3 -2 -1 0 1 2 3 4 5 ...23456
  // Depadded      -2 -1  0 1 2 3 4 5 6 ...01234 
  // 
  //                                T A ...CGATT 
  //                      N A G G G T A ...   mAlignStart = 0
  //                N  N  N A G G G T A ...   mAlignStart = -2
  //                          G G G T A ...   mAlignStart = 2
  if (PIdx < 0){
    //		cout<<"\tPIdx < 0"<<endl;
    return(PIdx+1);
  }else if (PIdx < mPadded2Depadded.size()){
    //cout<<"\tPIdx < "<<mPadded2Depadded.size()<<endl;
    // If it's in the range of consensus bases, neither before or after, 
    // then return it's position. 
    return(mPadded2Depadded[PIdx]);
  }else{
    //cout<<"\tPIdx >= "<<mPadded2Depadded.size()<<endl;
    // PIdx must be beyond the range of consensus bases.  So we simply 
    // add the difference.  In the example above, suppose we want to 
    // convert the padded coordinate of the T in CGATT, which is 16, to 
    // depadded.  The last depadded coordinate will be 12.  The diff 
    // in consensus size and PIdx will be 16-15+1=2.  The depadded idx 
    // of position 14 is 12.  12+2 = 14 as output.
    int lastPadded = mPadded2Depadded.size() -1;
    int diff = PIdx - lastPadded;
    return(mPadded2Depadded[lastPadded] + diff); 
  }  
}


//---------------------------------------------------------------------------
// * IsBase
//---------------------------------------------------------------------------
//
// Returns true of character is a valid base.
//
bool
Sequence::IsBase(char base)
{
  switch (base)
    {
      case 'A': 
      case 'C': 
      case 'G': 
      case 'T': 
      case 'a': 
      case 'c': 
      case 'g': 
      case 't': 
	return(true);
	break;
      default:
	return(false);
    }
}


//---------------------------------------------------------------------------
// * IsBaseNX
//---------------------------------------------------------------------------
// Returns true of character is a valid base, an n, or an x
//
bool
Sequence::IsBaseNX(char base)
{
  switch (base)
    {
      case 'A': 
      case 'C': 
      case 'G': 
      case 'T': 
      case 'a': 
      case 'c': 
      case 'g': 
      case 't': 
      case 'x':
      case 'X':
      case 'N':
      case 'n':
	return(true);
	break;
      default:
	return(false);
    }
}


//---------------------------------------------------------------------------
// * ComplementBase
//---------------------------------------------------------------------------
void 
Sequence::ComputeGCStats(int &NumBases, int &NumGC, int &NumCPG){
  NumBases = size();
  NumGC = 0;
  NumCPG = 0;
  ToLowerCase();
  int ssize = (*this).size();
  for(int i = 0;i< ssize;i++){
    if (((*this)[i] == 'c') || ((*this)[i] == 'g')) NumGC++;
    // Make sure out test doesn't go past end. 
    if (i <= ((*this).size()-2)){
      if (((*this)[i] == 'c') && ((*this)[i+1] == 'g')){
	NumCPG++;
      }
    }
  }
}


//-------------------------------------------
// * OnlyUpperCaseBases
//-------------------------------------------
// If s contains only elements of {a,c,g,t,n,x,A,C,G,T,N,X}, then this 
// function returns false if a character of s is anything except A,C,G, or T. 
//
// NOTE: Notice that this ONLY works if input is already a valid 
// sequence, not arbitrary ASCII.   
// 
bool
Sequence::OnlyUpperCaseBases()
{
  register int ssize = size();
  for(register int i = 0;i< ssize;i++){
    if ((*this)[i] & 0x28) return(false);
  } 
  return(true);
}

//-------------------------------------------
// * OnlyUpperCaseBases
//-------------------------------------------
// If s contains only elements of {a,c,g,t,n,x,A,C,G,T,N,X}, then this 
// function returns false if a character of s is anything except A,C,G, or T. 
//
// NOTE: Notice that this ONLY works if input is already a valid 
// sequence, not arbitrary ASCII.   
// 
bool
Sequence::OnlyUpperCaseBases(int start,int cklength)
{
  register int ckend = start+cklength;
  for(register int i = start;i< ckend;i++){
    if ((*this)[i] & 0x28) return(false);
  } 
  return(true);
}


//---------------------------------------------------------------------------
// * ContainsOnlyBasesSafe
//---------------------------------------------------------------------------
bool
Sequence::ContainsOnlyBasesSafe()
{
  for(register int i = 0;i< size();i++){
    if (!IsBase((*this)[i])) return(false);
  }
  return(true);
}


//---------------------------------------------------------------------------
// * ContainsOnlyBasesSafe
//---------------------------------------------------------------------------
bool
Sequence::ContainsOnlyBasesSafe(int start,int cklength)
{
  register int ckend = start+cklength;
  for(register int i = start;i< ckend;i++){
    if (!IsBase((*this)[i])) return(false);
  }
  return(true);
}


//---------------------------------------------------------------------------
// * ContainsNonBaseSafe
//---------------------------------------------------------------------------
bool
Sequence::ContainsNonBaseSafe(int start,int cklength)
{
  register int ckend = start+cklength; 
  for(register int i = start;i< ckend;i++){
    if (!IsBase((*this)[i])) return(true);
  }
  return(false);
}


//---------------------------------------------------------------------------
// * ComplementBase
//---------------------------------------------------------------------------
char
Sequence::ComplementBase(char base)
{
  switch (base){
    case 'A': return('T');
    case 'C': return('G');
    case 'G': return('C');
    case 'T': return('A');
    case 'a': return('t');
    case 'c': return('g');
    case 'g': return('c');
    case 't': return('a');
  }
  return(base);
}

//---------------------------------------------------------------------------
// * ToLowerCase
//---------------------------------------------------------------------------
// 
void
Sequence::ToLowerCase()
{
  int ssize = size();
  for(int i = 0;i<ssize;i++){
    (*this)[i] = tolower((*this)[i]);
  } 
}

//---------------------------------------------------------------------------
// * ToUpperCase
//---------------------------------------------------------------------------
// 
void
Sequence::ToUpperCase()
{
  for(int i = 0;i<size();i++){
    (*this)[i] = toupper((*this)[i]);
  } 
}

//---------------------------------------------------------------------------
// * ReverseComplement
//---------------------------------------------------------------------------
// 
void
Sequence::ReverseComplement()
{
  // Reverse the sequence
  reverse(begin(),end());
  register int ssize = size();
  for(register int i = 0;i<ssize;i++){
    (*this)[i] = ComplementBase((*this)[i]);
  } 
}

//---------------------------------------------------------------------------
// * Complement
//---------------------------------------------------------------------------
// 
void
Sequence::Complement()
{
  register int ssize = size();
  for(register int i = 0;i<ssize;i++){
    (*this)[i] = ComplementBase((*this)[i]);
  } 
}


//---------------------------------------------------------------------------
// * Nify
//---------------------------------------------------------------------------
// Converts anything in the sequence that isn't a base into an n. 
//
void
Sequence::Nify()
{
  register int ssize = (*this).size();
  for(int i = 0;i<ssize;i++){
    if(!IsBase((*this)[i])){
      (*this)[i] = 'n';
    }
  }
}

//---------------------------------------------------------------------------
// * Reverse
//---------------------------------------------------------------------------
// 
void
Sequence::Reverse()
{
  // Reverse the sequence
  reverse(begin(),end());
}


//---------------------------------------------------
// * Print
//---------------------------------------------------
//
void
Sequence::Print(ostream &os)
{
  os<<(*this).mName<<" "<<(*this).mDescription<<endl;
  for(int i = 0;i<size();i++){
    os<<(*this)[i];
  }
  os<<"\n";
}


//---------------------------------------------------
// * PrintSeq
//---------------------------------------------------
//
void
Sequence::PrintSeq(ostream &os)
{
  for(int i = 0;i<size();i++){
    os<<(*this)[i];
  }
  os<<"\n";
}


//---------------------------------------------------
// * Find
//---------------------------------------------------
// Returns the location of target within this sequence
// or -1 if target does not occur in this sequence. 
// This version accepts sequence iterators, allowing 
// us to search for some subset of the input string.  
//
long
Sequence::Find(Sequence::iterator iStart,Sequence::iterator iEnd)
{
  Sequence::iterator location;

  location = search(begin(),end(),iStart,iEnd);

  if (location == end()){
    return(-1);
  }else{
    long offset = location - begin();
    return(offset);
  }
}

//---------------------------------------------------------
// * FindEndOfTag
//---------------------------------------------------------
// 
//
void
Sequence::FindEndOfTag(Sequence &TargetSeq,
		       int TagStart,int TagEnd,
		       long &Location,bool &bIsReversed)
{
  bIsReversed = false;

  // Sanity check on input. 
  if (TargetSeq.end() < (TargetSeq.begin()+TagStart)) return;
  if (TargetSeq.end() < (TargetSeq.begin()+TagEnd)) return;

  // Ensure that both sequences are in lower case. 
  TargetSeq.ToLowerCase();  
  this->ToLowerCase();

  // Find the location of the of the tag. (Returns -1 for failed match)
  Location= this->Find(TargetSeq.begin()+TagStart,TargetSeq.begin()+TagEnd);
  
  // Location will be start of the tag.  We want the location
  // of the last character in the tag, so add the tag length. 
  if (Location >=0){
    Location = Location + (TagEnd-TagStart);
  }

  // If that failed to find a match, try matching the reverse compliment 
  // of the Target sequence. 
  if (Location < 0){
    TargetSeq.ReverseComplement();
    Location= this->Find(TargetSeq.end()-TagEnd-1,TargetSeq.end()-TagStart-1);
    // Location reflects the end of the tag sequence so no adjustment is 
    // needed
    bIsReversed = true;
    TargetSeq.ReverseComplement(); // Restore sequence to original sense.
  }
}

//---------------------------------------------------
// * Find
//---------------------------------------------------
// Returns the location of target within this sequence
// or -1 if target does not occur in this sequence. 
//
long
Sequence::Find(Sequence &target)
{
  Sequence::iterator location;

  location = search(begin(),end(),target.begin(),target.end());

  if (location == end()){
    return(-1);
  }else{
    long offset = location - begin();
    return(offset);
  }
}


//---------------------------------------------------------
// * FindEndOf
//---------------------------------------------------------
// Searches for an occurence of TargetSeq or it's reverse 
// compliment and returns the location of the end of that 
// sequence.   
//
long
Sequence::FindEndOf(Sequence TargetSeq, bool &bIsReversed)
{
  long location;

  // ensure that both target and source have same case. 
  TargetSeq.ToLowerCase();  
  ToLowerCase();

  location = Find(TargetSeq.begin(),TargetSeq.end());
  
  // Location will be start of the tag.  We want the location
  // of the last character in the tag, so add the tag length. 
  if (location >=0){
    location = location+TargetSeq.size()-1;
    bIsReversed = false;
  }

  // If that failed to find a valid set of endpoints, try matching the
  // reverse compliment of the FA sequence. 
  if (location<0){
    TargetSeq.ReverseComplement();
    location = Find(TargetSeq.begin(),TargetSeq.end());
    if (location >=0){
      location = location+TargetSeq.size()-1;
      bIsReversed = true;
    }
  }
  return(location);
}

//---------------------------------------------------
// * Print
//---------------------------------------------------
//
void
SequenceVector::Print(ostream &fout)
{
  for(int i = 0;i<(*this).size();i++){
    fout<<(*this)[i].mName<<" "<<(*this)[i].mDescription<<endl;
    (*this)[i].Print(fout);
  }
}


//---------------------------------------------------
// * FindEndOf
//---------------------------------------------------
// Finds the end of the query by searching all of the sequences
// in the vector.   Looks for both fwd and reverse compliment
// of the Query Sequence. 
//
void
SequenceVector::FindEndOf(Sequence &QuerySeq,int &SeqIdx,
			  long &SeqLocation,bool &bIsReversed)
{
  // Compare query to each sequence in this vector in turn. 
  long loc = -1;
  int sIdx = 0;
  int sEnd = (*this).size();

  do{
    loc = (*this)[sIdx].FindEndOf(QuerySeq,bIsReversed);
    sIdx++;
  }while((loc<0) && (sIdx < sEnd));

  SeqLocation = loc;
  SeqIdx = sIdx-1;  
}



//---------------------------------------------------
// * CopyToString
//---------------------------------------------------
void
Sequence::CopyToString(string &sout)
{
  register Sequence::iterator s = this->begin();
  register int i;
  register int Size = this->size();

  for(i = 0;i<Size;i++){	 
    sout+=*(s+i);
  }
}


//---------------------------------------------------
// * CopyToCString
//---------------------------------------------------
void
Sequence::CopyToCString(char *cstring)
{
  register Sequence::iterator s = this->begin();
  register int i;
  register int Size = this->size();

  for(i = 0;i<Size;i++){	 
    cstring[i]=*(s+i);
  }
  cstring[i]=0;

  //int i;
  //for(i = 0;i<this->size();i++){
  //  cstring[i]=(*this)[i];
  //}
  //cstring[i]=0;
}


//---------------------------------------------------
// * InitFromString
//---------------------------------------------------
// Assumes the string is an actual sequence. 
void
Sequence::InitFromString(string &source)
{
  StripWhiteSpace(source);
  for(int i = 0;i<source.length();i++){
    push_back(source[i]);
  }
}


//---------------------------------------------------
// * ConcatenateCString
//---------------------------------------------------
// Same as ConcatenateString but for c-language 
// strings (i.e. char*)
//
void
Sequence::ConcatenateCString(char *cstring)
{
  register int clen = strlen(cstring);
  register int oldSize = size();
  register long idx;

  // Make space for it.  
  resize(oldSize+clen);
  for(idx = 0;idx<clen;idx++){
    (*this)[oldSize+idx] = cstring[idx];
  }
}

//---------------------------------------------------
// * ConcatenateString
//---------------------------------------------------
// Assumes the s is an actual sequence.  This
// Is a pretty basic operation, performed many tens of 
// thousands of times by file reading functions, so it 
// should be as fast as possible. 
//
void
Sequence::ConcatenateString(string &s)
{
  register int slen = s.size();
  register int oldSize = size();
  register long idx;

  // Make space for it.  
  resize(oldSize+slen);
  for(idx = 0;idx<slen;idx++){
    (*this)[oldSize+idx] = s[idx];
  }
}


//---------------------------------------------------
// * ConcatenateNonWhiteSpace
//---------------------------------------------------
// Is a pretty basic operation, performed many tens of 
// thousands of times by file reading functions, so it 
// should be as fast as possible.  This variant exists
// to simplify filtering out whitespace, without some
// call to something like StripWhiteSpace, which is slow. 
//
void
Sequence::ConcatenateNonWhiteSpace(string &s)
{
  // Maximum size it might be. 
  register int slen = s.size();
  register long idx;
  register int newSize = size(); 
  register char base;
  
  // Make space for it.  
  resize(newSize+slen);
  for(idx = 0;idx<slen;idx++){
    base = s[idx];
    if ((base != ' ') && (base != '\n') && (base != '\t')){
      (*this)[newSize++] = base;			
    }
  }
  // newSize will equal the size of the concatenated sequence plus
  // one.  The Sequence has been resized to include all possible 
  // white space, so need to shrink it back down to the actual 
  // size, newSize -1. 
  resize(newSize-1);
}



//---------------------------------------------------
// * ConcatenateStringSafe
//---------------------------------------------------
// Checks each letter to see if it is a base before 
// adding. Is a pretty basic operation, performed many 
// tens of thousands of times by file reading functions, 
// so it should be as fast as possible. 
//
// 
void
Sequence::ConcatenateStringSafe(string &s)
{
  register int slen = s.size();
  register int oldSize = size();
  register long idx;
  register char c;
  // Make space for it.  
  resize(oldSize+slen);
  for(idx = 0;idx<slen;idx++){
    c = s[idx];
    if (IsBase(c)){
      (*this)[oldSize+idx] = c;
    }else{
      //cout<<"Non base where base expected."<<endl;
      //exit(1);
    } 
  }
}

//---------------------------------------------------
// * Concatenate
//---------------------------------------------------`
// Adds Seq to the end of this sequence.  Name and other
// info are unchanged. 
//
void
Sequence::Concatenate(Sequence &Seq)
{
  register int clen = Seq.size();
  register int oldSize = size();
  register long idx;

  // Make space for it.  
  resize(oldSize+clen);
  for(idx = 0;idx<clen;idx++){
    (*this)[oldSize+idx] = Seq[idx];
  }  
}

//---------------------------------------------------
// * InitFromCString
//---------------------------------------------------
// Assumes the cstring is an actual sequence.
void
Sequence::InitFromCString(const char *cstring)
{
  register int clen = strlen(cstring);
  register long idx;

  // Make space for it.  
  resize(clen);
  for(idx = 0;idx<clen;idx++){
    (*this)[idx] = cstring[idx];
  }
}

//---------------------------------------------------
// * InitFromCString
//---------------------------------------------------
// Assumes the cstring is an actual sequence.
void
Sequence::InitFromCString(char *cstring)
{
  register int clen = strlen(cstring);
  register long idx;

  // Make space for it.  
  resize(clen);
  for(idx = 0;idx<clen;idx++){
    (*this)[idx] = cstring[idx];
  }
}

/*

//---------------------------------------------------
// * GetSubSequence
//---------------------------------------------------
// Returns in Seq precisely the bases from idx start to 
// idx end, inclusive. 
// 
void
Sequence::GetSubSequence(Sequence &Seq,long start,long end)
{
  // Clear any previous sequence in the target.
  Seq.erase(Seq.begin(),Seq.end());

  if (end < 0) return;

  // Truncate ends
  if (start < 0) start = 0;
  if (end >= this->size()) end = this->size()-1;
  
  if ((end - start) >= 0){
	 for(long idx = start;idx <= end;idx++){
		Seq.push_back((*this)[idx]);
	 }

	 string description = (*this).mDescription;
	 // KJD Need to root these character arrays out. Suggest writing some
	 // string utilities that handle int->string, float->string, etc. 
	 char buf[500];
	 sprintf(buf,"[%ld,%ld]",start,end);
	 description = description+buf;
	 Seq.mDescription = description;
	 Seq.mName = (*this).mName;
  }
}
*/

//---------------------------------------------------
// * GetSubSequence
//---------------------------------------------------
// Returns in Seq precisely the bases from idx start to 
// idx end, inclusive. 
// 
void
Sequence::GetSubSequence(Sequence &Seq,long start,long end)
{
  // Clear any previous sequence in the target.
  Seq.erase(Seq.begin(),Seq.end());

  if (end < 0) return;

  // Truncate ends
  if (start < 0) start = 0;
  if (end >= this->size()) end = this->size()-1;
  
  long ssize = end - start + 1; 
  if (ssize > 0){
    Seq.resize(ssize);
    long subidx = 0;
    for(long idx = start;idx <= end;idx++){
      Seq[subidx++] = (*this)[idx];
    }

    string description = (*this).mDescription;
    // KJD Need to root these character arrays out. Suggest writing some
    // string utilities that handle int->string, float->string, etc. 
    char buf[500];
		// +1 to convert it to one based, since that is what the outside
		// world will want for coordinates. 
    sprintf(buf,"[%ld,%ld]",start+1,end+1);
    description = description+buf;
    Seq.mDescription = description;
    Seq.mName = (*this).mName;
  }
}

//---------------------------------------------------
// * Trim
//---------------------------------------------------
// 
//
void
Sequence::Trim(long tstart,long tend)
{
  for(int q = 0; q <= tstart;q++){
    (*this)[q]='X';
  }

  for(int q = tend; q < (*this).size();q++){
    (*this)[q]='X';
  }
}




//---------------------------------------------------
// * AddSequence
//---------------------------------------------------
// Copies the given sequence into the vector. 
//
void
SequenceVector::AddSequence(Sequence &Seq)
{
  /// KJD  Why the Hell would I make this function and 
  // not just do (*this).push_back(Seq) ?????


  // Increase the size of the vector to accomodate new sequence
  GrowVec();
  if (Seq.size() >0){
    register int ssize = Seq.size();
    for(register int idx = 0;idx < ssize;idx++){
      (*this)[size()-1].push_back(Seq[idx]);
    }
  }
  (*this)[size()-1].mDescription = Seq.mDescription;
  (*this)[size()-1].mName = Seq.mName;
}

//---------------------------------------------------
// * AddSubSequence
//---------------------------------------------------
// Copies the given subsequence into the vector. 
//
void
SequenceVector::AddSubSequence(Sequence &SourceSeq,long start,long end)
{
  // Increase the size of the vector to accomodate new sequence
  GrowVec();
  SourceSeq.GetSubSequence((*this)[size()-1],start,end);
  string description = SourceSeq.mDescription;
  // KJD Need to root these character arrays out. Suggest writing some
  // string utilities that handle int->string, float->string, etc. 
  char buf[500];
	// +1 to convert it to one based, which is what outside world
	// will want to see. 
  sprintf(buf,"[%ld,%ld]",start+1,end+1);
  description = description+buf;
  (*this)[size()-1].mDescription = description;
  (*this)[size()-1].mName = SourceSeq.mName;
}

//---------------------------------------------------
// * GrowVec
//---------------------------------------------------
//
void
SequenceVector::GrowVec()
{
  Sequence temp;
  temp.mDescription = "Empty";
  this->push_back(temp);
}


//---------------------------------------------------
// * TotalBases
//---------------------------------------------------
//
long
SequenceVector::TotalBases()
{
  long sum=0;
  for(int i = 0;i<this->size();i++){
    sum+=(*this)[i].size();
  }
  return(sum);
}


// I want reverse sort, so define my own compare
bool lt_myint(int c1,int c2){return(c2<c1);};

//---------------------------------------------------
// * N50
//---------------------------------------------------
// N50 is the max length L such that 50% of all bases 
// lie in contigs of size >= L
//
int 
SequenceVector::N50()
{
  vector<int> sizeList;
  
  GetSortedSizes(sizeList);
  
  int totalBases = TotalBases();
  int halfBases = (int)((float)totalBases /(float) 2.0);
  int sizeSum = 0;
  for(int i = 0;i< sizeList.size();i++){	
    sizeSum += sizeList[i];
    if (sizeSum >= halfBases) return(sizeList[i]);
  }

	// This should never happen. 
	return(-1);
}


//---------------------------------------------------
// * GetSortedsizes
//---------------------------------------------------
//
void 
SequenceVector::GetSortedSizes(vector<int> &sizeList)
{
  for(int i = 0;i< size();i++){
    sizeList.push_back((*this)[i].size());
  }
  sort(sizeList.begin(),sizeList.end(),lt_myint);
}

//---------------------------------------------------
// * MaxSize
//---------------------------------------------------
//
int
SequenceVector::MaxSize()
{
  int maxS = 0;
  for(int i = 0;i< size();i++){
    if ((*this)[i].size() > maxS) maxS = (*this)[i].size();
  }
  return(maxS);
}



//---------------------------------------------------
// * SizeStats
//---------------------------------------------------`
// 
void
SequenceVector::SizeStats(float &mean,float &variance)
{
  double SizeSum=0,DeviationSum = 0,n=(double)this->size();
  int seqNum;

  for(seqNum = 0; seqNum < this->size();seqNum++){
    SizeSum +=  (*this)[seqNum].size();
  }
  mean = SizeSum/n;

  for(seqNum = 0; seqNum < this->size();seqNum++){
    double deviation = pow((double)(*this)[seqNum].size()-mean,(double)2.0);
    DeviationSum += deviation;
  }
  variance = (1.0/(n-1.0))*DeviationSum;
}


//==========================================================================
// ------ | SequenceQuality|-------------------------------------
//==========================================================================

//---------------------------------------------------
// * Reverse
//---------------------------------------------------`
// 
void
SequenceQuality::Reverse()
{
  // Reverse the sequence
  reverse(begin(),end());
}


//---------------------------------------------------
// * Concatenate
//---------------------------------------------------`
// Adds Qual to the end of this quality.  Name and other
// info are unchanged. 
//
void
SequenceQuality::Concatenate(SequenceQuality &Qual)
{
  register int clen = Qual.size();
  register int oldSize = size();
  register long idx;

  // Make space for it.  
  resize(oldSize+clen);
  for(idx = 0;idx<clen;idx++){
    (*this)[oldSize+idx] = Qual[idx];
  }  
}

//---------------------------------------------------
// * Print
//---------------------------------------------------
//
void
SequenceQuality::Print(ostream &fout)
{
  int i;
  fout<<">"<<mName<<" "<<mDescription<<"\n";
  for(i = 0;i<size();i++){
    fout<<(*this)[i]<<" ";
    if ((!(i%50)) && (i!= 0)) fout<<"\n";
  }
  // Don't want to print out two new lines if by chance the 
  // last index happens to be a multiple of 50
  if ((i-1)%50) fout<<"\n";
}

//---------------------------------------------------
// * Print
//---------------------------------------------------
//
void
SequenceQuality::Print()
{
  Print(cout);
}

//---------------------------------------------------
// * Countphred20
//---------------------------------------------------`
// 
int
SequenceQuality::CountPhred20()
{
  int count = 0;
  for(int baseNum = 0;baseNum< this->size();baseNum++){
    if ((*this)[baseNum] >= 20){
      count++;
    }
  }
  return(count);
}

//---------------------------------------------------
// * FindPhred20
//---------------------------------------------------`
//
long 
SequenceQuality::FindPhred20(long start)
{
  for(long bNum = start; bNum < (*this).size();bNum++){
    if ((*this)[bNum] >= 20) return(bNum);
  }
  return((*this).size()-1);
}

//---------------------------------------------------
// * RFindPhred20
//---------------------------------------------------`
//
long 
SequenceQuality::RFindPhred20()
{
  return(RFindPhred20((*this).size()-1));
}


//---------------------------------------------------
// * RFindPhred20
//---------------------------------------------------`
//
long 
SequenceQuality::RFindPhred20(long start)
{
  for(long bNum = start; bNum >=0 ;bNum--){
    if ((*this)[bNum] >= 20) {
      return(bNum);
    }
  }
  return(0);
}

//---------------------------------------------------
// * GrowVec
//---------------------------------------------------`
//
void
SequenceQualityVector::GrowVec()
{
  SequenceQuality temp;
  this->push_back(temp);
}

//---------------------------------------------------
// * Trim
//---------------------------------------------------`
// 
//
void
SequenceQuality::Trim(long tstart,long tend)
{
  for(int q = 0; q <= tstart;q++){
    (*this)[q]=0;
  }

  for(int q = tend; q < (*this).size();q++){
    (*this)[q]=0;
  }

}

//---------------------------------------------------
// * Print
//---------------------------------------------------`
// 
//
void
SequenceQualityVector::Print(ostream &fout)
{
  for(int i = 0;i< size();i++){
    (*this)[i].Print(fout);
  }
}


//---------------------------------------------------
// * AddQuality
//---------------------------------------------------
// Copies the given sequence into the vector. 
//
void
SequenceQualityVector::AddQuality(SequenceQuality &Qual)
{
  // Increase the size of the vector to accomodate new sequence
  GrowVec();
  for(int idx = 0;idx < Qual.size();idx++){
    (*this)[size()-1].push_back(Qual[idx]);
  }
  (*this)[size()-1].mDescription = Qual.mDescription;
  (*this)[size()-1].mName = Qual.mName;
}

//---------------------------------------------------
// * Phred20Stats
//---------------------------------------------------`
// 
long
SequenceQualityVector::Phred20Total()
{
  long PhredSum=0;

  for(int seqNum = 0; seqNum < this->size();seqNum++){
    int phred20 = (*this)[seqNum].CountPhred20();
    PhredSum += phred20;
  }
  return(PhredSum);
}

//---------------------------------------------------
// * Phred20Stats
//---------------------------------------------------`
// 
void
SequenceQualityVector::Phred20Stats(float &mean,float &variance)
{
  double PhredSum=0,DeviationSum = 0,n=(double)this->size();
  vector<float>p20Vec;
  int seqNum;

  for(seqNum = 0; seqNum < this->size();seqNum++){
    float phred20 = (*this)[seqNum].CountPhred20();
    p20Vec.push_back(phred20);
    PhredSum += phred20;
  }
  mean = PhredSum/n;

  for(seqNum = 0; seqNum < this->size();seqNum++){
    double deviation = pow((double)(p20Vec[seqNum]-mean),(double)2.0);
    DeviationSum += deviation;
  }
  variance = (1.0/(n-1.0))*DeviationSum;
}


